Effects of image charges, interfacial charge discreteness, and surface roughness 
on the zeta potential of spherical electric double layers 
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We investigate the effects of image charges, interfacial charge discreteness, and surface roughness 
on spherical electric double layers in electrolyte solutions with divalent counter-ions in the setting 
of the primitive model. By using Monte Carlo simulations and the image charge method, the zeta 
potential profile and the integrated charge distribution function are computed for varying surface 
charge strengths and salt concentrations. Systematic comparisons were carried out between three 
distinct models for interfacial charges: 1) SURF1 with uniform surface charges, 2) SURF2 with 
discrete point charges on the interface, and 3) SURF3 with discrete interfacial charges and finite 
excluded volume. By comparing the integrated charge distribution function (ICDF) and potential 
profile, we argue that the potential at the distance of one ion diameter from the macroion surface is 
a suitable location to define the zeta potential. In SURF2 model, we find that image charge effects 
strongly enhance charge inversion for monovalent interfacial charges, and strongly suppress charge 
inversion for multivalent interfacial charges. For SURF3, the image charge effect becomes much 
smaller. Finally, with image charges in action, we find that excluded volumes (in SURF3) suppress 
charge inversion for monovalent interfacial charges and enhance charge inversion for multivalent 
interfacial charges. Overall, our results demonstrate that all these aspects, i.e., image charges, 
interfacial charge discreteness, their excluding volumes have significant impacts on the zeta potential, 
and thus the structure of electric double layers. 

PACS numbers: 82.70.Dd, 02.70.Uu, 61.20.Ja, 61.46.Df 



I. INTRODUCTION 

The study of charged interfaces in electrolytes is a 
problem of fundamental importance to biophysics, elec- 
trochemistry, and colloidal science^— A proper under- 
standing of the electric double layer (EDL) structure is 
essential to predict the stabilization of colloidal disper- 
sions and the properties of biological systems. Under ap- 
propriate physical and chemical conditions, charged in- 
terfaces display complex and counter-intuitive phenom- 
ena such as the charge inversion and like charge at- 
traction, which attracts a great theoretical and exper- 
imental interest £ These phenomena have been exten- 
sively observed in different systems including DNAs, self- 
assembled membranes and colloidal particles^ - — 

In the generally accepted (by chemists) Gouy- 
Chapman-Stern theory; 11 ' 12 an EDL is composed of an 
internal Stern layer, where some counter-ions are tightly 
bound to the charged interface, and an outer diffuse layer, 
where counter-ions exert thermal motions. The ion dis- 
tributions in the diffuse layer are usually calculated us- 
ing the Poisson-Boltzmann (PB) theor y 13 ! 14 Being of 
mean field nature, PB ignores the excluded-volume ef- 
fects as well as electrostatic correlations of ions. It is 
popularly believed that PB fails in the presence of mul- 
tivalent ions or highly charged interfaces! 15 ' 16 Various 
methods, including but not limited to, modified PB the- 
ories, integral equation theories, and density functional 
theory have been developed to describe physics beyond 
the Poisson-Boltzmann framework. For example, it is no- 



table that the state-of-the-art classical density functional 
theory^— has incorporated accurate the hard-sphere re- 
pulsion and electrostatic correlation up to high precision, 
yielding results that are in quantitative agreement with 
Monte Carlo simulations. 

While most of previous studies on the EDL model 
interfacial charges as uniformly distributed over the 
macroion surfaces, in reality, the surface charges are 
better modeled as discrete. When colloidal particles 
were placed in electrolyte solutions, the surface chemical 
groups release hydrogen cations to the solvent, resulting 
in negatively charged interfacial ions, whose strength de- 
pends on the environment conditions. As an extremal 
example, phospholipids in aqueous solutions can carry 
a variable charge between — 4e and +le under different 
physiological conditions Recently, there has been 

many works demonstrating 2 ^— that interfacial charge 
discreteness has an important influence on the microion 
distribution near flat interfaces. Faraudo et aL 28 ' 29 dis- 
covered that discrete interfacial charges can lead to an 
inversion of selectivity observed in a protein channel in 
the presence of multivalent cations. Fleck et al<^ showed 
that a (on average) negatively charged membrane consist- 
ing of mobile surface groups can adsorb a like-charged 
DNA. These motivate a systematical study on the effect 
of discrete interfacial charges. 

Charged objects immersed in electrolytes usually have 
much lower dielectric constants than water. Therefore 
polarization charges (image charges) on the interfaces are 
a relevant issue in the study of EDL. The understand- 
ing of image effects is of recent interests for electrostatic 
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interactions of soft matter systems ! 22 ' 31 " — Counter-ions 
that are attracted to the interface are repelled by their 
likely charged image charges when they approach the in- 
terface. Therefore image charges reduce the ion density 
in the vicinity of charged interface. This effect however 
diminishes as the surface charge strength increases, as 
shown by Torrie et alM- long time ago. For strongly 
charged surface, image charges push the whole EDL out- 
ward by a small distance, but otherwise has no significant 
influence on the phenomena such as charge inversion, see 
recent works by Wang and Ma. 22 ' 33 ' 34 For recent reviews 
on image charge effects in colloidal and biological sys- 
tems, see also references! 37 ' 38 We note, however, there is 
a recent work by Boda et al^ showing that image charge 
effects may play important roles in biological systems. 

If the interfacial charges are discrete, they also have 
image charges. In aqueous solvents and for planar in- 
terfaces, these images are almost identical to the source 
charges and therefore essentially double the surface 
charge density. This may lead to substantial modification 
of EDL structure, as we shall show in the present work. 
Finally, charged interfaces are usually not smooth at mi- 
croscopic scales. Surface roughness at atomic scales may 
interfere with image charge effects, and therefore change 
the physics of charge inversion phenomena. It is our pur- 
pose to study the interplay between image charges, inter- 
facial charge discreteness, and surface roughness using a 
model system of spherical colloid. 

To study image charge effects for a generic interface, 
it is necessary to numerically solve the corresponding 
boundary value problem for Poisson equation. This is 
usually too time-consuming to be feasible in Monte Carlo 
simulations. For some of recent works on smooth inter- 
faces, see references ! 40 ' 41 In this work, we shall explore 
the effects of image charges, interfacial charge discrete- 
ness, and surface roughness on the zeta potential and the 
charge inversion phenomenon for spherical geometries by 
using a recent method of multiple imagesi 42 ' 43 We shall 
compare three different toy models for charged interfaces: 
a smooth surface with uniform surface charge density, 
a smooth surface with discrete interfacial charges, and 
a rough surface with discrete interfacial charges. This 
problem is difficult to address using conventional method 
of spherical harmonic expansion ; 44 ' 45 due to the intensive 
computation cost of the polarization potential. This dif- 
ficulty can be surmounted using the recently developed 
method of multiple images i 42 ' 43 

The term zeta potentia l '' 47 is intimately related to the 
Smoluchowski theories for electrophoresis: It is defined 
as the electrical potential in the interfacial double layer 
(DL) at the location of the supposed slipping plane versus 
a point in the bulk fluid away from the interface. The 
existence of a slipping plane is one of the fundamental 
assumptions of electrophoresis theories, but we are not 
aware of any direct experimental evidence for its exis- 
tence. Defined as such, zeta potential can not be directly 
measured, but can only be inferred indirectly from elec- 
trokinctic data, through the application of Smoluchowski 



theories. Numerical calculations for the zeta potential 
have been made under different kind of geometries and 
various electrolytes ! 48 ' 49 For 1:1 electrolyte with low salt 
concentration, it was shown that the mean-field PB the- 
ory has a very good agreement with the simulation re- 
sults. But in other cases, for example, with a 2 : 2 salt 
in electrolytes, the PB fails to predict the result both 
qualitatively and quantitatively^ 

For strongly charged surfaces, there is a layer of 
counter-ions strongly bound to the dielectric interface. 
If inversion of electrophoresis mobility occur, this con- 
densed layer must move with the colloid in electrophore- 
sis. Diehl and Levin 5 - 1 - argued that in numerical computa- 
tion of zeta potential using Monte Carlo simulation, the 
slipping plane should be identified at about one counter- 
ion diameter away from the colloidal surface. In the 
presence of charge inversion, an alternative, but prob- 
ably even more natural choice of slipping plane would 
be the peak of the integrated charge distribution func- 
tion (ICDF). This later choice would rigorously identify 
the charge inversion with the reversal of colloidal mobil- 
ity, which generally assumed to be true. It is interesting 
to note that in most numerical simulations, the peak of 
ICDF is indeed about one micro-ion diameter away from 
the surface. We note further that near the charge in- 
version threshold, the precise location of slipping plane 
is only of minor importance, because the potential pro- 
file changes very slowly near the peak of ICDF. We shall 
follow the choice of Diehl and Levin 5 -^ in this work. 

The remaining of this work is organized as follows. In 
Section II, we present the three distinct models for in- 
terfacial charges and the simulation method used in this 
work. In Section III, we present the simulation results 
and discuss in details the effects of interfacial charge dis- 
creteness, image charges and surface roughness. Finally, 
we draw the concluding remarks in Section IV. 



II. MODELS AND METHODS 

We consider a charged colloidal particle with radius 
a = 2nm and dielectric constant £; = 2, hereafter re- 
ferred as the macroion, immersed in a 2 : 2 symmetric 
electrolyte. Such a macroion can be used to model mi- 
celles, dendrimers and other colloids i 5 - 2 - - — The aqueous 
solvent is modeled as a dielectric continua with a dielec- 
tric permittivity e = 78.3, while the ions are modeled 
as small hard spheres with diameter r = OAnm and with 
all charges in their centers. The spherical Wigner-Seitz 
(WS) cell model 5 - 7 - is employed for the boundary in our 
simulation. The macroion is located at the center of the 
cell (with radius R) and has a bare charge Qm = Zm^i 
surrounded by the solvent and microions. The electrolyte 
is treated at the level of restricted primitive model. The 
microions are confined in the spherical WS cell. There 
are N + counter-cations with valence Z + = 2 and 7V_ co- 
anions with valence Z_ = — 2 in the system. The whole 
system is charge-neutral, hence N+Z++N-Z-+Zm = 0. 
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FIG. 1: Schematics for three models of interfacial charges studied in this work: (a) SURF1: uniform surface charge. It is 
equivalent to a single point charge at the center of macroion; (b)SURF2: discrete interfacial point charges on the surface; (c) 
SURF3: discrete interfacial charges with excluded volume. Also in (a), we show the image charges of a mobile microion qi as 
blue empty circles. In (b) and (c), a discrete interfacial charge go produces image charges. The Kelvin image coincides with 
the source charge itself, while other images lie inside the sphere. 



To explore the effects of image charges, interfacial 
charge discreteness, and surface roughness on the struc- 
ture of EDL, we shall systematically compare three differ- 
ent models for interfacial charges, which arc illustrated in 
Fig. [TJ In model SURF1, the interfacial charges are uni- 
formly distributed on the colloidal surface. Equivalently, 
we can also put all surface charges onto the center of the 
spherical colloid. Because of spherical symmetry, the im- 
age charges of all interfacial charges cancel each other, so 
they have no influence on mobile microions in the bulk. 
In model SURF2, there are No point charges of valence 
Zd distributed on the smooth spherical colloidal sur- 
face. The total bare charge of the macroion is therefore 
Qm = QdNd- Three different valences of discrete inter- 
facial charges are studied: Zjj E { — 1,-2,-3}. In prac- 
tice, Zd = — 1 corresponds to the carboxylatc functional 
group^; Zd = —2 can be some kind of phospholipid 
such as DMPAj£2 finally the so-called PIP2 lipids^^ 
can have a valence between —2 and —4 depending on the 
solvent environment. The locations of these Njj charges 
are determined by running a MC minimization of the 
electrostatic energy near zero temperature with the con- 
straint that the charges remain on the colloidal surface. 
Typically, after 10 5 Nd MC moves, the interfacial charges 
form a Wigncr crystal on the surface. Model SURF3 
is the same as model SURF2, with the only difference 
that each interfacial charge acquires an excluded vol- 
umes (bumps) with diameter tjj = 0.4nm. For simplicity 
of computations, these bumps are assumed to have the 
same dielectric constant as the solvent. These two mod- 
els have also been used in the study of ion binding to 
polyclcctroly tcs 

The bumps on the interface in SURF3 may have two 
competing impacts. Firstly they increase the minimal 
distance between interfacial charges and mobile ions, by 
an amount of 7t>/2 = 0.2nm. This tends to weaken 



the interaction energy between interfacial charges and 
the mobile ions, and therefore suppress binding between 
them. Secondly the bumps also provide more interfa- 
cial area, so that more than one counter-ions can be 
adsorbed near a given interfacial charge. Competition 
between these two tendencies plays an important role in 
the local structure of EDL. 

The electric potential distribution. 3>(r), for a snapshot 
ion distribution satisfies the Poisson equation (Gauss 
unit), 

- V-e(r)V$(r) =4vr^<5(r-r 3 ), (1) 

i 

subject to standard electrostatic boundary conditions. 
The dielectric function e(r) takes £i inside the sphere 
and e D outside, S(-) is the Dirac delta function, and the 
index j runs over both the mobile microions in the bulk 
as well as the interfacial charges on the macroion sur- 
face. Numerically solving this Poisson equation in three 
dimensional space is time costly. Luckily for spherical ge- 
ometry, there is an efficient image charge algorithm.] 42 ! 43 
Given the dielectric boundary, the electrical potential 
at r generated by one unit point charge at r' is given by 
the electrostatic Green's function G(r, r') that satisfies 
the following differential equation, 

_ J_ V- £ (r)VG(r,r')=5(r-r'), (2) 

47T£ 

with homogeneous electrostatic boundary conditions. We 
use light italic font r' = |r'| to represent the radial dis- 
tance in spherical coordinates. For interfacial charges 
r' = a while for mobile ions r' > a. The Green's function 
is a linear superposition of the point-charge Coulomb po- 
tential in free space, 
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and the potential of all image charges due to the polar- 
ization of the macroion G- lm (r, r') ( with a unit source 
charge fixed at r'): 



G(r,r') = Go(r,r') + G im (r,r'). 



(4) 



For spherical geometry, the image potential was studied 
in detail in reference |42| . Here we invoke the result di- 
rectly: 



G im (r,r') = 



- r ya/r' 
■ - r K \ 



-fv(x/r K ) v 1 /a 



dx, (5) 



where v = e j {e\ + e Q ) and 7 = {e\ - e a )/(si + e Q ) w -1. 
The first term is due to a (likely charged) point im- 
age with charge —^ajr' at tk = r'a 2 /r' 2 (Kelvin im- 
age), while the second term is due to an (oppositely 
charged) line image extending from vk back to the center 
of sphere. These image charges are overall neutral. By 
contrast, a planar interface only produces a point image 
(Kelvin image), and is not charge neutral. The impor- 
tance of the line image is controlled by the ratio between 
the Debye length and the sphere radius, A/a. For our 
system, this parameter is not big, neither is it negligibly 
small, see Table Q] 

Using the /-point Gauss-Legendre quadrature to ap- 
proximate the line integral, Gi m (r, r') can be rewritten 
as the potential due to a total of / + 1 image charges^ 
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{w m , s m , 1 < m < 1} are the /-point Gauss weights and 
locations on the interval [—1,1]. Note that m = corre- 
sponds to the Kelvin image charge, i.e., we have ujq = —2 
and so = — 1. In this simulation we choose 7 = 2. The 
image charges for mobile and interfacial charges in mod- 
els are illustrated in Fig. 1 (a) and (b)(c), respectively. 

The total Hamiltonian of the system can be expressed 
as a sum of three contributions^ 
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(7) 



The first part £/j Mm is the interaction between macroion 
and microions, while the second part U™ m is the interac- 
tion between microions. The third part [/hs is the hard 
sphere repulsions between mobile ions, interfacial ions, 
the macroion, and the WS cell shell. It takes the positive 
infinity when any volume exclusion constraint is violated 
and zero otherwise. 

For SURF1, the interfacial charges are uniformly dis- 
tributed on the sphere, hence their image charges, when 
added up, cancel each other, and have no interaction with 
the mobile ions. Therefore C/j Mln is the direct Coulomb 
interaction between the central charge at the origin and 
the mobile ion, 
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(8) 



where Ib = e 2 /4ireoE kBT is the Bjerrum length, eo is the 
vacuum permittivity, and fcgT is the thermal energy, /3 = 
l/k B T. For SURF2 and SURF3, U^ m is the interaction 
between interfacial ions and the mobile ions, 
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(9) 



where G is given by Eq. ((4]), while r„ is the position of 
the nth interfacial ion. It should be noted that due to the 
symmetry of the Green's function, G(r i; r n ) = G(r„,ri), 
the interactions between the interfacial ions and the im- 
age charges of the mobile ions have been already in- 
cluded. Similarly, the microion-microion interaction (the 
in Eq. ([7])) can be written as, 
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(10) 

where Sij is the Kronecker delta. When i = j, it repre- 
sents the interaction of the charge and its images. 

Canonical-ensemble Monte Carlo simulations based on 
the standard Metropolis acceptantce/rejection rul o 64 i 65 
are carried out to obtain the equilibrium properties of 
the model systems. The initial configuration of each sys- 
tem is generated by randomly placing the ions into the 
simulation cell satisfying the constraints of the hard-core 
repulsion. The total number of mobile ions N varies 
from about 200 to 300, depending on the salt concen- 
trations and the WS cell radius. In each simulation, 
we perform 1.2 x 10 6 iV MC moves per particle. The 
first 10 5 ./V MC moves are performed using the tempering 
techniqu e) 66 ' 67 , where we start with T = 2100if (Bjerrum 
length Ib = O.lnm), and then slowly cool down to room 
temperature T = 298A" (Bjerrum length Ib = 0.71nm). 
This is followed by another 10 5 iV MC moves for equili- 
bration. Finally, 10 6 iV MC moves are performed to store 
the data for statistical analysis. In order to achieve better 
sampling efficiency, the acceptance ratio was kept about 
0.3 by adjusting the maximum step size of ion motion.— 
We find that the autocorrelation function (ACF) of to- 
tal energy typically decays to zero after about 12000 MC 
steps (for the whole system). Since in our simulation, 
about 200 x 10 6 MC steps are calculated, we have nearly 
20000 independent samples. The ACF for one typical 
run is illustrated in Fig. [2] 

We chose symmetric and multivalence 2 : 2 electrolytes 
instead of more common asymmetric 2 : 1 electrolytes 
mainly from the consideration of computational capac- 
ity. Nevertheless, 2 : 2 divalent salt systems are also 
experimentally relevan t 69 ' 70 and have been studied both 
analytically 7 ^— and numerically . 74 ' 75 It is known that 
the structure of EDL is mainly determined by valence 
of counterionsi 48 ' 49 In our simulations, the salt concen- 
trations of the systems take three values: G = 0.1, 0.2 
and 0.5M; the surface charge density of the macroion a 
varies from —0.1 to —0.6 Cm~ 2 , correspondingly the bare 
charges range from — 30e to — 180e, the Gouy-Chapman 
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FIG. 2: Energy autocorrelation function (ACF) as a function 
of MC steps. System parameters: the surface charge density: 
Q.lC/m 2 , Zd = 3, salt concentration C = Q.1M. For other 
parameters, see Table [I] 



length (Iqc = 1/2ttZ + Ib<j) ranges from 0.031 to 0.19nm. 
The radius of the WS cell R is chosen to be at least 16 
times longer than the Debye length A, so that the influ- 
ence of the cell wall on the EDL structure is negligible. 
The radius of the spherical macroion is kept a = 2nm. 
All simulation parameters are summarized in Table |U 

TABLE I: Relevant system parameters used in the Monte 
Carlo simulations 



As discussed in the Introduction, the zeta potential is 
defined as the average potential one micro-ion diameter 
away from the colloidal surface, £ = <J>(a + r). For spher- 
ical EDLs, this potential can be explicitly obtained by 
integrating the Poisson equation^ 




where Pi(r) is the mean density of the ith ion species, 
and Zi its valence. The integrated charge distribution 
function (ICDF) Q(r) as a function of radius r is given 
by 

Q{r) = Q M + [Z+N + (a,r) + Z_N-(a,r)), (12) 

where N±(a,r) arc the average numbers of posi- 
tive/negative ions in the spherical shell between a and 
r, and Qm the bare charge of macroion. It is found that 
Q(r) changes sign at r « z + 3r/4 when charge inversion 
takes place. The maximum of the ICDF curve is defined 
as the inverted charge, which equals zero if no charge 
inversion happens. See Fig. 01 Fig. [5] and Fig. [5] for illus- 
trations of ICDF. We have also calculated the standard 
deviation of zeta potential to make sure that the simu- 
lation results are accurate. For the system parameters 
shown in Fig. [2] for example, the zeta potential is found 
to be ( = 31.833mV, while its standard deviation is: 
0.279ml/. 



III. RESULTS AND DISCUSSION 
A. SURF1: Uniform surface charge distribution 



Si = 2 Colloidal dielectric constant 

e = 78.3 Solvent dielectric constant 

Zm = —30 ~ —180 Macroion valence 

Z± — ±2 Counterion and coion valence 

Zo = —1,-2,-3 Interfacial ion valence 

a = 2nm Macroion radius 

t — OAnm Microion diameter 

td = OAnm Diameter of interfacial ions in SURF3 

T = 298K Room temperature 

Ib = 0.71nm Bjerrum length 

Igc = 0.19 ~ 0.031nm Gouy-Chapman length 

C = 0.1,0.2,0.5M Three salt concentrations 

A = 0.48, 0.33, 0.21?im Corresponding Debye lengths 

R — 7.72, 6.13, 4.19n?n Corresponding WS cell radii 



The zeta potential and the integrated charge distri- 
bution function (ICDF) are calculated for three differ- 
ent surface charge models, SURF1, SURF2 and SURF3. 



In model SURF1, interfacial charges are smoothly dis- 
tributed on the spherical colloidal surface. The image 
charge potential of interfacial charges is averaged out and 
has no influence, but microions do have image charges 
that are inside the colloidal sphere. These images repel 
the source charge away from the colloidal surface, re- 
sulting in a depletion zone for mobile microions near the 
interface. It has been shown that except this depletion 
effect, image charges of the mobile ions have no influ- 
ence on the charge inversion phenomenal In agreement 
with these results, our simulation results show that the 
image charges only slightly change the value of the zeta 
potential. Illustrated in the first row of Fig. [3] are the 
zeta potential for model SURF1 with and without image 
charges, for three different concentrations of 2:2 salt. 

For comparison we also show in the bottom row of 
Fig. [3] the average potential $(a + r/2) at distance of 
one microion radius t/2 away from the interface (which 
is often used in literature as zeta potential) . We find that 
this potential is drastically different from $(a + r). More 
importantly, $(a + r/2) remains negative for the whole 
range surface charge strength studied for the cases of salt 
concentration 0.1M and 0.2M, and therefore shows no 
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FIG. 3: Top row: The zeta potential (average potential at one micro-ion diameter away from the interface) for SURF1 as the 
function of surface charge density a for both with and without image charge effects. Three different concentrations of 2:2 salt 
are calculated: (a) C — O.lAf, (b) C — 0.2M, and (c) C = 0.5A1". Bottom row: average potential at half micro-ion diameter 
away from the interface. 
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FIG. 4: Model SURF1. The integrate charge distribution 
functions (ICDFs). a = -0.6C/m 2 and C = 0.2M. The 
peak of these curves are at slightly more than one micro-ion 
diameter outside the surface. 



sign of charge inversion. By strong contrast, the potential 
$(a + t) indicates charge inversion at a ps — 0.35Cm -2 
(salt 0.1M) and a ps -0.3Cto" 2 (salt 0.2M). The 
plot for integrated charge distribution function Fig. |4] 
also clearly shows charge inversion for C = 0.2, <j = 
— 0.6C/m 2 . We therefore conclude that <£>(a + r), rather 
than $(a + r/2), is a good definition of the zeta poten- 
tial in term of the characterization of charge inversion 
phenomenon. This agrees with the results by Diehl and 
Levin^i 



B. SURF2: Discrete surface charges without 
volume effect 



In model SURF2, intcrfacial charges are point charges 
distributed on the sphere. The corresponding zeta poten- 
tial as a function of surface charge density is illustrated 
in Fig. [5j Three cases for the valences of interfacial 
charges are shown, with Zd = —1,-2 and —3, respec- 
tively. For Zd = — 1 (Fig.[5j left column), the zeta poten- 
tials with image charge effect become substantially lower 
than those with no image charge effect. Comparing with 
Fig. [31 it is evident that the image charge effects greatly 
enhance the overcharging tendency. This is clearly due 
to the effective doubling of interfacial charge valence by 
their image charges. 

Interestingly, in the cases of higher valences of interfa- 
cial ions, Zd = —2 (Fig. El center column) and Zd = —3 
(Fig. 03 right column), image charge effects influence zeta 
potential profiles in the opposite direction. The zeta po- 
tentials with the image charges remain uniformly below 
those without images, therefore the tendency of charge 
inversion is substantially suppressed by the image charge 
effects. Furthermore, the image charge effects grow with 
increasing surface charge density, which is the opposite 
of what was discovered for uniformly charged surfaces 
by Torrie et al.S^ It is also interesting to note that the 
suppression of charge inversion is more significant for 
Zd = — 2 than for Zd = — 3. 

These effects of image charge effects is also evident 
in the ICDF plots. In Fig. [5] we show ICDF curves 
for various interfacial charge valences, with fixed surface 
charge density a = — 0.6C/to 2 and salt concentration 
C = 0.2M. It is clear that for Z D = -1, the inte- 
grated charges with image charges always lies above that 
without image charges, while for Zd = —2 and —3, the 
curves with image charges lie below those without im- 
age charges. In fact, for these latter two cases, the image 
charge effects completely eliminate charge inversion. The 
suppression is again more pronounced for Zd = — 2 than 
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FIG. 5: The zeta potential for SURF2 as a function of mean surface charge density, with various valences of interfacial charges, 
both with and without image effects. Left: Zd = — 1; middle: Zd = —2; right: Zd = —3. It is clear from these plots that 
image charges enhance charge inversion for Zd = — 1, and suppress charge inversion for Zd = —2, —3. 
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FIG. 6: Model SURF2. The integrate charge distribution 
functions (ICDFs) for different valences of interfacial ions. 
a = -0.6C/m 2 and C = 0.2M. 



for Z D = -3. 

The huge suppression of charge inversion for divalent 
interfacial charges is likely due to the strong binding be- 
tween these interfacial charges and the divalent counter- 
ions, which completely neutralizes the interfacial charge 
group. Other counter-ions can no longer be adsorbed 
onto the nearby area if such a binding occurs. The bind- 



ing energy is about Zd£b/(t/2) « 7ksT without image 
charges and is further doubled by the image charge ef- 
fects. This binding may strongly decreases the effective 
surface charge density and therefore suppresses charge in- 
version. To verify this physical picture, we have looked at 
typical simulation snapshots of ion distributions near the 
colloidal surface, for model SURF2 (no bumps, with im- 
age charges). One typical simulation snapshot is shown 
in Fig. [7] for each of three case Zd = — 1, —2, —3. It is 
found that for Zd = — 2 (Fig. [3 center panel), about 90% 
of interfacial charges are closely bound to a counter-ions, 
and become completely neutralized. Such a binding is 
clearly strengthened by image charge effects. 

By contrast, for Zd = — 1 (Fig. [3 left panel), about 
46% of interfacial charges are closely bound to a counter- 
ions. Since the counter- ions carry charge 2e while in- 
terfacial charges carry — e, each of these bindings con- 
tribute to the charge inversion. Indeed model SURF2 
with Zd = — 1 goes charge reversal at a much lower sur- 
face charge density, comparing with SURF1, as one can 
see from Fig. [5] and Fig. 02 Finally, for Z D = -3 (Fig. [3 
right panel), we can see that the majority of interfacial 
charges are also bound to counter-ions. In this case, how- 
ever, one-to-one binding does not invert the charge, nor 
does it completely neutralize the interfacial charge group. 
Some interfacial charges actually can bind to two coun- 
terions. Their numbers are however not enough to in- 
vert the charge of the whole colloid. Model SURF2 with 
Zd = — 3 goes charge reversal at higher surface charge 
density than model SURF1. 

We note that the binding between interfacial charges 
and counterions are important for Zd = —2, —3 even 
in the absence of image charges. Image charge effects 
however substantially enhance these bindings. Our sim- 



FIG. 7: Typical snapshots of model SURF2 after equilibrium is achieved. Gray spheres: counterion; yellow spheres: coion; 
small white spheres: interfacial charges. Left: Zd — — 1; Center: Zd = —2; Right: Zd = —3. System parameters: salt 
concentration 0.2M, surface charge density 0.6C/m 2 . 



ulation demonstrates the following: comparing with uni- 
form surface charge distribution, interfacial charge dis- 
creteness enhances charge inversion if interfacial charge 
groups carries have smaller valence than counterions, and 
suppresses charge inversion if interfacial charge groups 
have equal or larger valence than counterions. Image 
charges strengthen these effects by effective doubling of 
surface charge density. 



C. SURF3: Discrete surface charges with finite 
exclusion volumes 

In model SURF3, there is a spherical excluded vol- 
ume (bump) with diameter td = OAnm around each 
interfacial charge, which prevents the counter-ions from 
getting closer than r = OAnm to the interfacial charge. 
Zeta potentials for SURF3 for various values of C and 
Zjj are shown in Fig. [8j Comparing with Fig. [5j we see 
that the effects of image charges become much smaller, 
in fact almost negligible, in the presence of volume exclu- 
sion effects of bumps. Nevertheless, for all three cases of 
interfacial charge valences, image charge effects slightly 
enhance charge inversion, in strong contrast with model 
SURF2. Furthermore, while in SURF2, the charge inver- 
sion threshold of surface charge density increases with 
the valence of interfacial charges, in SURF3, the charge 
inversion threshold actually decreases with the valence 
of interfacial charges. Another striking effect is that, the 
more realistic model SURF3 gives a much higher zeta po- 
tential than model SURF1, as one can see by comparing 
Fig. [5] with Fig. [3J These surprising results are proba- 
bly due to the extra interfacial area around the bumps, 
where multiple counter-ions can be attracted by a given 
surface charge. This effect becomes more important as 
the valence of interfacial charge increases. 

We also plot the ICDF curves in Fig. |H]for fixed salt 
concentration C = 0.2M and the surface charge density 
a = — 0.6 C/m . The figure shows the charge inversion is 
moderately enhanced by the image charges by examining 
the strength of the inverted charges. The peaks of ICDF 



curves are also pushed a little further away from the sur- 
face by the image charges, in agreement with previous 
studies by other groups i 22 ' 33 ' 34 ' 36 ' 37 

We also compare the zeta potential of models SURF2 
and SURF3 (both with image charges) in Fig. [TO] We 
find that for Zd = — 1, the data for SURF3 arc uniformly 
and substantially lower than those for SURF2, indicat- 
ing that charge inversion is strongly suppressed by the 
excluded volumes. By strong contrast, for Zd = —2 or 
—3, the zeta potential for SURF3 are substantially higher 
than those for SURF2, indicating that charge inversion 
is strongly enhanced by the bumps. 

What is the physics underlying these huge differences 
between SURF2 and SURF3? Recall that in model 
SURF3, each surface charge group has an exclusion vol- 
ume ("bump") with radius 2 A. These bumps have two 
competing effects: Firstly the minimal distance between 
surface ions and the ions in the bulk is A A in SURF3, 
instead of 2A in SURF2. Hence the maximal interaction 
energy between surface ions and counterions are approx- 
imately reduced by half in SURF3. This seems to be the 
dominant factor in the case of Zd = — 1 , leading to strong 
suppression of charge inversion in SURF3. Secondly the 
bumps provide extra adsorption area for the counterions. 
If the interaction between interfacial charges and coun- 
terions are already strong enough, these extra areas will 
lead to more counterions condensed near the interface, 
and therefore enhance charge inversion. This seems to 
be the dominant factor for the cases of Zd = —2 and 
Z D = -3. 



IV. CONCLUDING REMARKS 

We have studied the effects of image charges, inter- 
facial charge discreteness, and surface roughness on the 
zeta potential profile and the integrated charge distri- 
bution function for a strongly charged spherical colloid. 
Such a study is possible owing to the recent development 
of image charge methods for a spherical boundary, which 
approximates the image potential with a few point image 
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FIG. 8: The zeta potential for SURF3 as a function of surface charge density in cases both with and without image charges 
using two different definitions of the shear plane location. The valence of interfacial ions Zd = — 1, —2, —3 from left to right. 
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FIG. 9: Model SURF3. The ICDF curves for Z D = -1,-2 
and —3 for the cases with and without image charges. Salt 
concentration C = 0.2M and surface charge density a — 
-0.6C/m 2 . 



the nature of surface charge distribution. For uni- 
form surface charge (SURF1), the influence of im- 
age charges on zeta potential is minor, in agreement 
with previous results by other groups. For discrete 
surface charge groups with no excluded volumes 
(SURF2), we find that image charges strongly en- 
hance charge inversion if the surface charge groups 
are monovalent, and strongly suppress charge inver- 
sion if the surface charge groups are multi- valence. 
For discrete surface charge groups with excluded 
volumes (SURF3), we find that the effects of image 
charges only slightly enhance charge inversion. 



(3) Model SURF3 with discrete interfacial charges and 
finite excluded volumes gives a much higher zeta 
potential than model SURF1, where the interfa- 
cial charges are continuous and the colloidal sur- 
face is smooth. Finally, the effects of excluded vol- 
umes (bumps) of interfacial charges are substan- 
tial and depend crucially on the valence of surface 
charges. For monovalent surface charges, we find 
that surface roughness strongly suppress charge 
inversion, while for multi-valence surface charges, 
surface roughness strongly enhance charge inver- 
sion. These results show that roughness of charged 
interface is an extremely important factor. 



charges. Our main results are summarized as below: 

(1) In agreement with Diehl and Levin^i we find that 
the zeta potential defined at about one counter-ion 
diameter away from the colloidal interface provides 
a good indicator for charge inversion phenomenon. 

(2) The effects of image charges depend crucially on 



Overall, our simulation results demonstrate intricate 
and competing effects associated with image charges, in- 
terfacial charge discreteness, and surface roughness. All 
these factors can substantially affect the profile of the 
zeta potential and charge inversion. Therefore short scale 
details of charged interfaces need to be better clarified 
before the structure of EDL can be understood properly. 
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FIG. 10: The zeta potential for SURF2 and SURF3 (both with image charge effects) as a function of mean surface charge 
density. The valence of interfacial ions Zd = — 1, —2, —3 from left to right. 
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